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We study the diffusion of monochromatic classical waves in a disordered acoustic medium by 
scattering theory. In order to avoid artifacts associated with mathematical point scatterers, we 
model the randomness by small but finite insertions. We derive expressions for the configuration- 
averaged energy flux, energy density, and intensity for one, two and three dimensional (ID, 2D 
and 3D) systems with an embedded monochromatic source using the ladder approximation to the 
Bethe-Salpeter equation. We study the transition from ballistic to diffusive wave propagation and 
obtain results for the frequency-dependence of the medium properties such as mean free path and 
diffusion coefficient as a function of the scattering parameters. We discover characteristic differences 
of the diffusion in 2D as compared to the conventional 3D case, such as an explicit dependence of 
the energy flux on the mean free path and quite different expressions for the effective transport 
velocity. 
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I. INTRODUCTION 



The ongoing interest in the field of classical waves in 
complex media is caused by the importance of detection 
and imaging techniques that are based on wave prop- 
agation and scattering. This ranges from electromag- 
netic waves in optical and near infrared tomography jl| 
and microwave radars Id] to acoustic waves in ultrason- 
ics 13 and geophysics Complexity is often associated 
with inhomogeneities that cause scattering which com- 
plicates most imaging processes considerably. However, 
when used cleverly, the scattered field can also be used 
to improve imaging j^. Although length scales (with re- 
spect to the wavelength) and the degree of the disorder 
may vary considerably from field to field, methods and 
results have been shown to be interchangeable without 
much difficulty ^6]. Recent topics of interest include lo- 
calization of classical waves 13, 0, the transition from 
ballistic to diffusive wave propagation acoustic time- 
reversal imaging p^ . etc. Direct simulation by the exact 
solution of a well-known Helmholtz wave equation for a 
given realization of the medium is often the method of 
choice for given applications. The drawbacks of the brute 
force computational approach are the limited system size 
and statistics that can be achieved with given computer 
resources as well as the difficulty to distill general princi- 
ples out from the plethora of output data. The need for 
simple models with transparent results therefore remains. 

An analytic theory of wave propagation in disordered 
media necessarily relies on simple model scatterers, for 
which point scatterers, i.e. (re gula rized) (5-functions in 
real space, are often chosen jllllijl- Unfortunately, the 
scattering response of a single point scatterer can be- 
come non-causal, a pathological behavior that can not 
be solved by a simple momentum cutoff Especially 
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for the study of the frequency dependence over a wider 
range it is therefore necessary to use more realistic model 
scatterers. 

In this paper we wish to study a simple but not unre- 
alistic experiment for the determination of the scattering 
properties of scalar waves in a disordered bulk material. 
A signal is emitted by a source and detected by a receiver, 
both embedded in the medium at sufficiently large dis- 
tances from the boundaries. Ultimately, we are interested 
in the detector signal caused by a pulsed (broadband) sig- 
nal emitted by the source. After a first arrival we then ex- 
pect the so-called coda that arrives at later times due to 
multiple scattering at the random scatterers How- 
ever, combining both the effects of multiple scattering 
and the full frequency dependence of the scattering pro- 
cesses renders an analytical treatment difficult without 
additional approximations, such as a complete neglect of 
the frequency dependence of the scattering amplitudes 
|l5j |. In order to understand how to justify certain ap- 
proximations and eventually find better ones, we have 
carried out a study of the frequency dependence of the 
scattering properties of random media. We concentrate 
on the steady state in the presence of strictly monochro- 
matic sources, which distinguishes the present work from 
related studies of the propagation of narrow band pulses 
|l6j | . As main results we obtain the frequency dependence 
of macroscopic effective medium properties like the mean 
free path and the diffusion constant that depend on the 
microscopic parameters of the random scatterers. 

When the ratio between source-receiver distance and 
mean free path is small, wave propagation is predom- 
inantly ballistic. When this ratio is large, energy and 
intensity propagation is governed by the diffusion equa- 
tion hi hj. Both these regimes are well understood. 
However, many imaging applications operate on length 
scales where the mean free path and the source-receiver 
distance are comparable. This is especially the case in 
geophysics where mean free paths range from a few hun- 
dred meters up to tens of kilometers [13 ■ The behavior at 
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this crossover regime between balhstic and diffuse wave 
(intensity) propagation is of considerable interest Q and 
also subject of the present study. 

Here we present an analytical formalism on monochro- 
matic wave intensity and energy propagation in one di- 
mensional, two dimensional and three dimensional (ID, 
2D and 3D) homogeneously disordered media using re- 
alistic model scatterers. We determine the relative con- 
tributions of diffusively and coherently propagated waves 
as the source-receiver distance increases. We did no find 
many theoretical studies of wave propagation in two di- 
mensional random media in the literature 

U m, al- 

though several experiments on quasi 2D systems have 
been carried out 0) ^| ■ Another possible test for our 
2D theory is comparison with numerical studies, which 
for very large systems are much cheaper than in the 3D 
case. 

The remainder of this paper is organized as follows. 
In sections II-IV we start with defining our model sys- 
tem and the basic equations, addressing the scattering 
matrices of single scatterers and discussing the average 
amplitude propagators in the frequency domain. The in- 
tensity, energy flux and energy density are discussed in 
section V. Results on the frequency dependence of the dif- 
fusion constant and its dependence on the model param- 
eters are discussed in section VI. In general, the results 
for ID systems are easily obtained, whereas our results 
for 3D systems agree with findings previously reported 
by others. The mathematics in the 2D case is not triv- 
ial, however, and the derivations are summarized in the 
appendix. We end with the conclusions. 

II. DEFINITIONS AND BASIC EQUATIONS 

A. Microscopic equations 

We describe the propagation of (scalar) acoustic waves 
in a microscopic model system. Specifically, we consider a 
ID, 2D or 3D acoustic medium with wave velocity cq and 
a mass density po- The medium contains n randomly dis- 
tributed scatterers per unit length, area or volume and 
we treat the dilute limit in which the average distance 
between scatterers is much larger than their radius a. 
The internal wave velocity of a scatterer is Cint and, for 
simplicity, the difference in mass density with the sur- 
rounding medium is disregarded. The waves are emitted 
by a monochromatic point source oscillating at frequency 
uj positioned at the origin. The wave amplitude ip(^: re- 
lated to the hydrostatic pressure by ~ dti^^j and to 
local particle velocity by = —p'^^Vtjjuj, then obeys the 
wave equation: 

(V2 + (r) d'i) (r; t) = -QpoS (r) cos {cut) . (1) 

The source term chosen here corresponds to a volume 
injection term. The source emits plane waves for ID, 
cylindrical waves for 2D and spherical waves for 3D me- 
dia. In all cases Q is in units of length per unit time. The 



wave velocity profile of the entire medium c (r) contains 
the information of the positions of the scatterers (in ID 
r = x). 

The Green function of the Hclmholtz equation in 
the real space and frequency domain reads 

[V2 + >,l-V (r; Lo)] G (r, r'; uj) = S {r - r') , (2) 

where kq = i-^/cq, the length of the wave vector in the ho- 
mogeneous medium. V (r; uj) is the scattering or impurity 
potential, a sum over all individual scattering potentials: 

N 

V{r;u;)^>,l{l-j-')J2Q{a-\r-r,\). (3) 

Q is the Heaviside step function, with 8 (x) =0 when 
X < and 1 otherwise. The velocity contrast is defined 
as 7 = Cint/ Co so that the single scatterer potential is 
"attractive" when 7 < 1 and "repulsive" when 7 > 1. 
Eq. Q describes a spherical potential, but the precise 
shape is not relevant when the scatterers are sufficiently 
small compared to the wave length. 

The amplitude of the wave field is related to the Green 
function as 

i;^ (r; t) - -Qpo Re {e-^"*G (r, r' = 0; c^)} . (4) 

The intensity (r; t) is the square of this expression. 
Related physical properties are the energy flux: 

F^{r;t) = -—dt4>^{r;t)\/^^{r;t), (5) 
Po 

and the energy density: 

(r; t)^^ (r; t)f + c'^ (r) (S^Vc. (r; t)f) , 

(6) 

recognized as the sum of the potential and kinetic energy 
contributions, respectively. For a monochromatic source 
with frequency u these observables contain a time inde- 
pendent contribution and a second term oscillating with 
frequency 2u;. We concentrate on the constant part by 
time-averaging over one period. Expressed in terms of 
the Green function this yields 

/..(r)-^|G(r,r' = 0;a;)^ (7) 



Fc. (r) = - Im {G (r, r' = 0; u:) VG* (r, r' = 0; c.)} , 

(8) 

T^.(r)-^(|VG(r,r' = 0;c.)|^ 

+ ^|G(r,r' = 0;^)|^). (9) 
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B. Macroscopic equations 

The properties of the wave field depend via the Green 
function on the exact configuration of scatterers. How- 
ever, in large systems, different realizations of the ensem- 
ble give similar responses (ergodicity) . The similarities 
in the response can be studied by calculating the configu- 
rational average. This average is the connection between 
the microscopic description and the macroscopic (effec- 
tive) medium properties. 

The macroscopic (diffusively scattered) intensity of 
pulsed sources is usually described by the diflFusion 
equation: 



dt{I{v-t))^D^^ {I{r;t)), 



(10) 



where the brackets denote the configuration average and 
D is the diffusion constant. In spite of neglecting the fre- 
quency dependence of the scattering processes, this ap- 
proximation is known to work well in certain cases (under 
the condition that the source receiver distance is much 
larger than the mean free path) ^E^. 

In order to obtain the steady-state diffuse intensity of a 
monochromatic wave field, the diffusion equation H1U|) is 
not sufficient. The energy density ( and not the intensity) 
of the wave field is the conserved property. Eq. H10() is 
therefore only valid if the intensity is strictly proportional 
to the energy density. In general, the averaged energy 
transport is governed by Pick's first law: 



(F^ [v)) = ~D{uo)V{W^(v)), 



ill) 



accounting for the frequency dependence of the diffusion 
constant. In the steady-state problem and outside the 
monochromatic source the proper Laplace equation is 



{W^ (r)) = 0. 



(12) 



In ID, the s-wave scattering condition corresponds to 
equivalence of for either refiection or transmission. 
R{lo) is the reflection coefficient at a step discontinuity, 
it can be obtained by imposing flux conservation across 
the scatterer boundary. This gives 



R{u) 



^i?o(l-e*'"»4a/7) 
1 - i?ge*''o4°/T 



(15) 



where i?o = (7 — 1) / (7 + !)■ In the same way we can 
derive an expression for the scattering matrix element of 
the s-wave channel 5*0 (related to the scattering phase 
shift 5o by Sq = exp (i2(5o)). In 2D 



So {uj) 



7 Jo (K0Q/7) H{ (noa) ~ Ji (Koa/7) Hq (koo) 
7J0 (Koa/7) -^1^^ (i^oa) - Ji (f^aa/l) J^o^'' (^oa) 



(16) 



In 3D the Bessel (Ji) and Hankel {h\^'^) functions are 

replaced by the spherical Bessel {ji) and Hankel {h^P) 
functions. The scattering matrix element then simplifies 
to I2II 



o ( X -»2Koa cot (Koa/7) -1-^7 

cot (Koa/7) — 17 



(17) 



IV. 



THE CONFIGURATION-AVERAGED 
PROPAGATOR 



Now we switch to the case of multiple scattering at the 
proposed model scatterers. The wave propagator in a dis- 
ordered medium after configuration averaging is dressed 
with a self-energy E. In reciprocal space it reads ^3 



III. SCATTERING MATRICES 

Here we discuss the properties of a single model scat- 
terer in the system {N = 1 in equation ||2J)). The response 
of a system containing a monochromatic source (in the 
origin), a receiver (at r) and a single "s-wave" scatterer 
(at Ti) can be expressed in terms of Green functions of 
the homogeneous system {V = 0) po| : 

G(r,r' = 0;c^) = Go (r;c^) 

+ Goi\r~r,\;Lu)toiiu)Go{r;uj). (13) 

This expression is valid in the far field limit (r, 3> A) 
and when scattering is isotropic (A ^ a), where A is the 
wavelength. 

The transition (t-) matrix elements for s-wave scatter- 
ing are related to the scattering matrix elements by 



to (t^) 



2iKo_R (uj) 

2i [So {lo) - 1) 

Itxik^^ {So {oj) 



(ID) 
(2D) 
(3D) 



(14) 



(G(k,k';c.)) = 



k2 -E(k;w) 



(27r)^ (k-k'), 



(18) 

where d is the dimension and S^'^' the Dirac delta func- 
tion. When n, the density of scatterers, is low, interfer- 
ence between multiply scattered waves by different sites 
may be disregarded. In this "single site approximation" 
the self-energy does not depend on k and it is simply 
given by 0| 



S {lu) — nto {lu) . 



(19) 



This approximation does not restrict the scattering 
strength since to is the full scattering matrix of the sin- 
gle scatterer. Interference effects from multiple scatter- 
ing at different scatterers cause localization that is known 
to be important in ID (where the localization length is 
of the order of the mean free path) and in 2D media 
(where the localization length is a transcendental func- 
tion of the mean free path). In 3D, localization can be 
disregarded except for very strong scattering media 0. 
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Here we restrict ourselves to purely non-localized trans- 
port phenomena, remembering that we can always find a 
region where this type of transport is dominant. 

Fourier transforming (|18|l with self-energy given by 
(|19|1 to real space gives the averaged Green function that 
depends only on the source-receiver distance (G (r; to) — 
{G{r,r' = 0;u!))). In ID the amplitude propagators are 
exponentially damped plane waves: 



Gi\x\;Lo) 



1 



^inc {uj)\x\ 



(20) 



in 2D they arc cylindrical: 



and in 3D spherical 



if ^>o , 
{-Ke (w) r) if w < 



-1 



(22) 



@. In Eq. (Enna Kg is the "renormalized" effective 
wave vector: 



sgn (uj) Kr (uj) + i 



21 f (^)- 



(23) 



Kr(yi) = |RcKe(cL;)| and t^^^ iuj) = 2|ImKe([j)|, the 
mean free path. We retrieve the Green functions for the 
homogeneous systems (Go) by letting n or to go to zero. 
Properties of the averaged response to a pulsed signal 
can be studied by calculating the Fourier transform to 
the time domain, as was done in Refs. and p^. 



V. THE CONFIGURATION-AVERAGED 
INTENSITY END ENERGY 

We derive here the configuration averaged intensity, 
energy flux and energy density in the frequency domain. 



A. The Bethe-Salpeter equation 



This is the Bethe-Salpeter equation in position space, 
where Ho is the coherent intensity (Ho = |(G)|^) and 
r is the irreducible vertex function. The lowest order 
approximation that still accounts for multiple scattering 
is 

r(ri,r2,r3,r4;w) =nr (cj) -J^'^' (ri - rg) 

x,5(^)(ri-r2)(5(^'(r3-r4). (26) 

This reduces the Bethe-Salpeter equation to 

n(r;c^) - Ho {t-uj) 

+ nr(c^) y"Aino(|r-ri|;c^)n(ri;c^). (27) 

In reciprocal space this integral equation becomes a geo- 
metric series that can be summed as 



n(/fc;tj) 



Ho ik-uj) 



l-nP (w)no (A:;w)' 



(28) 



In order to be able to calculate the Fourier transform 
of n (/c; w) , an expression for IIq (fc; w) is needed. It is cal- 
culated as the Fourier transform of the coherent intensity 
and this results in ID in 



Ho {h,uj) 



21) 



in 2D in 



({2ti,iff + l) ({klff-l 



(29) 



no(fc;w) = 



„2 arcsm , , 



= , (30) 



and in 3D in [12 



n,ik;u)='^ ^''Xr'\ (31) 
47r ktf 

The calculation of the vertex function F is discussed in 
the next subsection. 



Ensemble averaging the intensity {Tj) gives us 

(/^(r)) = ^n(r;c.), (24) 

where n(r;(jj) = {^|G(r,r' — Q\uo)\^^ is the average of 
the squared Green function propagator. It is given by 

n (r; u)) = Ho (r; ^) + j d'^r^d'^rd'^rsd'^ri (G (r, ri;uj)) 

X (G* (r, r2;u;)) T (n, r2, rs, r4; lj) 

X (G (rg, r' = 0; u) G* {n, r' = 0; u)) . (25) 



B. Energy conservation and the Ward identity 

It is well known that for a given approximation for the 
self-energy, one is not free to choose the vertex correc- 
tion. Here we take advantage of the flux conservation 
constraint to obtain F without additional microscopic 
calculations. The energy flux from the monochromatic 
source (on average) points outwards. In the steady state 
case the following condition must hold for the averaged 
flux in direction n: 

(n-F^(r)) cx ^n-?, (32) 
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where r is the unit vector in the radial direction. In ID 
this condition reads 

(F„ (x)) oc sgn (x) . (33) 

The microscopic expression for the average energy flux is 



(n-F^ (r)) 



X Im {{G (r, r' = 0; cj) n • VG* (r, r' ^ 0; lu))} 



Im n" (r;u;) , 



(34) 



which defines the function H". The vertex function is 

the same as for the intensity, so we can express H" in 
reciprocal space as 



n" (k;u;) 



(k;c.) 



1 -nr (w)no (fc;w) 



(35) 



Ilg {k; cu) is the coherent energy flux in direction n that 
is given by the Fourier transform of 



In ID it is straightforward to show that condition H33|l 
can only be fulfilled when Eq. H41|l is fulfilled as well. The 
Ward identities are relations between self-energy and ver- 
tex corrections. We can identify Eq. 1)41(1 as the Ward 
identity for our problem. We now have all the ingredi- 
ents to calculate the Fourier transform of H28|) and H35|) 
to calculate the averaged intensity and energy fiux re- 
spectively. 



C. Flux 

Using the Taylor expansions in the limit fc — s- 0, we 
find an expression for C (from Eq. H37|) ') in 2D and 3D: 



C 



dr Im {G (r; lu) drG* (r; lj)} r'^d'^ 
ino-i(fc = 0;c.)92no(fc;a;)|,^Q 



(42) 



The average fiux in ID is obtained by directly Fourier 
transforming (|35|l : 



(r; uj) = G (r; uo) n • VG* (r; u) . 



(36) 



In 2D and 3D the averaged microscopic expression for 
the energy flux should match the macroscopic condition 



Im|n"(r;c.)| = --g^n.?, 
which in reciprocal space reads 



(37) 



Re|n"(k;cj)| = - (n-k) 2"*-^^ (38) 

where G is real and depends on frequency and the model 
parameters. Hq (fc;w) is an even function of fc. We know 

how IIq (k; w) depends on fc, as the Fourier transform in 
2D reads 



nf^ (k;w) = - (n-k) 27 



X j drJi(kr)rG{r]uo)drG* {r;uo) , (39) 





and in 3D 

n" (k;t^) = - (n-k)47ri 



X j drji{kr)r^G{r-uj)drG* {r-uj). (40) 



The Taylor series of 11" (k; w) around fc = only contains 
odd terms. So, in the limit that fc ^ 0, condition H38() 
can only be fulfilled by Eq. (|35|) when 



nV (w) =no^(fc = 0;w). 



(41) 



1/ (2£;)' 



:sgn (x) 



(43) 



We show how to calculate G in 2D case in the ap- 
pendix. With the result, the projection of the average 
fiux becomes 



(n • (r)) = 



Q^Po arctan (2Kr^/) 
87r2 r 



n • r, 



(44) 



while in 3D it is straightforward to calculate G from Eq. 
(I42|l and the projection of the average fiux then reads 



(n • F^ (r)) 



{Any r2 



-n • r. 



(45) 



Letting if — > oo (k^ — * l^ol) recovers the fiux of a 
monochromatic source in an unperturbed medium. 

It is interesting to see that, in contrast to the 3D case, 
in ID and 2D the average flux depends on both the mean 
free path and the real part of the effective wave vector. 
So the scattering mean free path limits the energy flux 
in ID and 2D, but not in 3D. In a strongly scattering 2D 
medium, in which the wave energy is not (yet) localized, 
the dependence on the arctangent should be observable. 



D. Intensity 

The total average intensity is proportional to the prop- 
agator n [r]Lo)^ that can be obtained by calculating the 
Fourier transform: 



d'^k 



(46) 
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Ho {k) is given by ^ in ID, ^ in 2D and ^ in the 
3D case. In ID and 2D this integral diverges because 
in the steady state case with a monochromatic source, 
energy does not escape fast enough to infinity due to 
the scatterers. This is analogous to the fact that the 
Poisson equation (the diffusion equation in steady state 
with source term) for a line or planar source has no well- 
defined solution. 

The gradient of the intensity exists in all cases. In ID 
it is constant and the derivative of 11 {x; uj) is given by 



dxTl (x; Lo) — —sgn (x) 



1 



1 



(47) 



4£/ «2 + i/(2^^)^ 
The gradient of 11 in 2D is expressed as an integral by 



_„ / N ^ f dk 
Vn(r;[j) = -r / r- 



PJi [kr] 



{k-L0)-Ii^^{k = {)-Uj) 



rf (r;w) , 



(48) 



which defines a function / (r; lo), that represents the gra- 
dient in the r direction. We split it up into a coherent 
(co/i) and a "totally diffusive" {td) part and a crossover 
correction (cr) 

/ (r; u) = fcoh {r; oj) + ftd (r; uj) + U {r; uj) . (49) 

The coherent part is connected to the unscattered inten- 
sity, therefore 

fcoh ir]uj) = dr\G{r-uj)\'^ 

= -^Rc\^{Kr + i/2£f)H[^'^ {{Kr + i/2lf)r) 



X H^^'{{Kr-i/2lf)r)} 



In the appendix it is shown that 

„ / N arctan (2Kr^f ) 1 i « ^ 
ftd{r;uj) = __L_Z11Z_5-1(2k,^^), (51) 



TT'^2Krif r' 



with 

g{2Kref) = 1- 



{2Kr£f) 2Krif eLTCtan{2Kr£f)' 



(52) 



This part decays like 1/r, much slower than the coher- 
ent and crossover contributions. It is the part that de- 
scribes the intensity gradient when energy transport is 
completely governed by Pick's first law, that is why we re- 
fer to this term as the "totally diffusive" part. When the 
total gradient is approximated by the just the sum of the 
coherent and the totally diffusive contribution, the gra- 
dient first decays exponentially until the source-receiver 
distance is approximately two to three mean free paths 
and then the 1/r decay is dominant. However, in this 
approximation it is neglected that close to the source the 
diffusive field will look different than far away from the 



source. The third term of/, the crossover term, describes 
this difference. In the appendix it is shown that 



ftd {r;uj) + fcr {r;uj) 

OO 

= -J ^Jiikr) fc'n,e(fc;c.), 



where 



IIsc (k-Lu) 



Uo^{k = 0;uj)llo{k;uj) 
n„'{k;Lo)-n^'{k^O;L. 



(53) 



(54) 



We did not find an analytical expression for the integral 
(|53|l and thus we need to evaluate it numerically. The 
crossover term vanishes for r/£f or r/£f ^ 1 and 
peaks at r/£f « 0.3. Only around this value of r/£f the 
gradient (in absolute value) is overestimated significantly 
(up to 25%) when we approximate it by just the sum of 
coherent and "totally diffusive" terms. 

In 3D the Pourier transform (|46|l converges and the 
intensity is well-defined. We rewrite 

(55) 

where 



h{T/£f)^ d^ 



4(e + i) 



(2 (e+1)- In (1 + 2/0)' + 



- 1 



X e 



(56) 



We were not able to solve (jBfill analytically as well. In 
the 3D case, the intensity is a function of £y only (it 
does not depend on n^). Eq. (|55|l consists of three terms 
(n = IIco/i + Hid + T\cr)- The first term is proportional 
to the coherent intensity, the second term is the alge- 
braically decaying diffuse term (the only term that is not 
exponentially decaying). We plot the sum of the first and 
second term multiplied by r (in units of 1/ (IGtt^^/)) in 
figure J^l as a function of r/^/ (on the right axis). The 
third term is again the crossover correction to the total 
intensity if we approximate 11 by only the first two terms. 
A plot of the crossover correction divided by the sum of 
the first two terms is shown in figure ^ (left axis). The 
crossover term vanishes ioi r/£f or r/£f ^ 1 and 
peaks at r/£f « 0.3. It can thus be concluded that the in- 
tensity can very well be approximated by just the sum of 
coherently and totally diffusively propagated intensities, 
as the total intensity in 3D will never be overestimated 
by more than 5% using this approximation. 

To complete this discussion we show the final results 
for the gradient of the average intensity in ID and 2D: 



9. (/c (x)) 



sgn [x] 



(57) 



V (r)) « ?^ (/,„, (r; oj) + ftd {r; co)) . (58) 
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FIG. 1: r (n^oh + litd) in units of (IBtt^^/)"^ (dotted line, 
right axis) and Her / (Rcoh + ^td) (solid line, left axis) as a 
function of r/£f (the source-receiver distance in number of 
mean free paths) for the 3D disordered medium. 



where fcoh and ftd are given by (|5n|l and l(5T|l respectively. 
The average intensity in 3D is approximated well by 

We obtain these expressions from our first-principles cal- 
culations that enable us to study not only the ballistic 
and diffusive limits, but also the crossover regime when 
r/if « 1. From this we observe that we can approximate 
the average intensity well by just the coherent and diffu- 
sive contributions. Furthermore, we saw that already at 
r/£f « 0.3 the diffusive intensity is higher than the co- 
herent intensity. This does not mean that when a pulsed 
source is used, we should see signs of the crossover to 
the diffusive regime at this point, because the diffuse 
peak is much broader than the coherent peak, so this 
crossover point is at larger values of r/£f, as was previ- 
ously reported Obviously, our present model system 
has been assumed to be boundless. In a finite slab ge- 
ometry boundary scattering, which is beyond the scope 
of this study, would of course affect the results. 

E. Energy density 

To derive a first-principles expression for the diffusion 
constant from Fick's Law Hll|) . we still have to calculate 
the average energy density given by 

(W^.(r)) = ^((|VG(r,r' = 0;c.)|') 

+ (oj^c-^r)\G{r,r' = 0;ojf)). (60) 

The first term is the average potential energy density and 
the second term corresponds to the kinetic energy. 



We start with the potential energy term in 2D and 3D. 
We define 

(|VG(r,r' = 0;c.)|')=n(r;c.). (61) 

" 2 

The Fourier transform of no(= |(VG)r) diverges, which 
means that we can not use the same procedure as we 
used for the intensity and the flux. According to the 
Bethe-Salpeter equation 

n(r;c^) =no (r;c^)+n^i(fc = 0;u;) 

X y"dVno(ri;c^)n(|r-ri|;t^). (62) 

This integral diverges as well because of the strong sin- 
// 

gularities in IIq (also when the gradient is calculated in 
the 2D case). When averaging, scatterers are effectively 
moved around the medium and for every configuration 
the contribution to the total average response is calcu- 
lated. However, because of the stronger singularities in 
// 

Ho (remember that every scatterer becomes a new source 
of spherical waves) this is not possible when the receiver 
position coincides with a scatterer position. The rea- 
son for this is the point receiver assumption and the far 
field scattering approximation. We can circumvent this 
problem by omitting a small volume/area around ri with 
radius of approximately one wavelength. This slightly 
modifies the probability distribution function form "com- 
pletely random" to "non-overlapping" (with the receiver) 

// 

in order to avoid the divergencies. We then find that 11 
is given by: 

n(r;w) = |«:e|'n(r;w). (63) 

In principle, our original expression for 11 should now be 
multiplied by a factor exp {—ro/£f), where ro is the radius 
of omission, so as long as the mean free path is longer 
then a few wavelengths omitting this small volume does 
not influence the results. Furthermore, even if scattering 
is strong and the mean free path is of the order of the 
wavelength, this factor will not be of importance. 

The second term of Eq. H6()|l. the kinetic energy, can 
be split: 

{lu^c-' (r)|G(r,r' = 0;u;)|') 

= K^U (r; u;) - (v (r; c.) |G (r, r' = 0; c.) |') . (64) 

Now the condition that the scatterer position can not co- 
incide with the receiver position ensures that the second 
term vanishes, due to the step function in the potential 
0. We can thus just disregard this term. 
In ID it is straightforward to prove that 

^(|x|;u;) = |/^e|'^(|x|;c^), (65) 



8 



always holds. We have to impose the condition the the 
receiver can not coincide with a scatterer to make sure 
that 

(lo^c-^{x)\G{x,x' ^Q-Lo)^) ^ nln[\x\-uj). (66) 

Only under the restrictions mentioned here, the averaged 
energy density in ID, 2D and 3D can be expressed as 
being proportional to the intensity: 

{W^{v)) = ^^{\KA^ + ni){I^{v)). (67) 

and this thus means that only the gradient of the energy 
density is well-defined in the ID and 2D cases. 



VI. THE DIFFUSION CONSTANT 

Using the Bethe-Salpeter equation with the Ward iden- 
tity we find expressions for the average energy flux 
I45II . the (gradient of) the average intensity H57I59|I . The 
average energy density is just proportional to the aver- 
age intensity, Eq. ((HH). When r/^/ ^ 1 we expect (fTT|> 
to hold and as the gradient of the average energy density 
and the average flux are now known we find an expression 
for the diffusion constant from This means that the 
diffusion constant can be written as 



■D(t^) ^ ^Ceff {uj)£f (w) , 



where in the ID and 3D case 



eff (w) = Co 



2Kr \ko\ 



kI -f 1/ + 



(68) 



(69) 



and in the 2D case 



Co 



+ 1/ [21 fY + 



:g {2njf) , (70) 



where g {2Kr£f) is given by H52I) . The effective transport 
velocity in 2D reduces to (|69|l in the weak scattering limit. 

We can now investigate the frequency dependence of 
the diffusion constant for a medium with monodisperse 
scatterers. We relate the scatterer density n to the aver- 
age distance between scatterers ((ds)) so that n = {ds)~^ 
in ID, n = Att'^ (4)"^ in 2D and n = 3 (47r)~^ (4)"^ 
in 3D. Let us focus on the diffusion constant of the 2D 
medium. We write 



{ds 



toiaKo), (71) 



so that the dimensionless property OKg depends on the 
dimensionless frequency kqu (= Loa/co) and two dimen- 
sionless model parameters, i.e. the velocity contrast 
7(= Cint/co) and the average distance between scatter- 
ers in number of scatterer radii {{ds) /a). The real and 



T— - ■ 1 r— 




■ \ V ■ 


y=0.2 










y=0.9 










y=i.i 










y=10.0 


















1 





FIG. 2: Diffusion constant of the 2D disordered medium in 
units of Co a, as a function of the dimensionless frequency Koa 
for four different scatterer- medium velocity ratios (7). The 
scatterer density is determined by setting {ds) /a = IQ. 



imaginary parts of OKg are needed to obtain the diffusion 
constant 



(oKo) = |Re{aKe (aKo)}| , 



if (oKo) 



2 |Im {flKe (fl'*o)}| 



(72) 



(73) 



The diffusion constant for a 2D medium is plotted in 
Fig. 121 The relevant frequency range is from KQa{= 
uja/co)= to Koa « 7r/2, as for higher frequencies the 
isotropic scatterer assumption is no longer valid. For the 
plot, the density of scatterers was determined by setting 
{ds) I a = 10, increasing this value shifts the curves up. 
The shape of the curves is predominantly determined by 
the mean free path. The effective transport velocity Ce// 
only deviates considerably from co when the scatterer 
velocity and the frequency are small and the scatterer 
density high. For the diffusion constants shown in the 
plot, this is only the case when 7 = 0.2. This is also the 
only case that shows resonances in the relevant frequency 
range. Lowering the internal velocity of the scatterers 
even more, would "pull in" more resonances in the rele- 
vant frequency range. These resonances show up because 
of resonances in the mean free path. When the scatterer- 
medium velocity ratio is increased, the mean free path 
(and thus the diffusion constant) increases until the ra- 
tio is larger than unity and it will drop again. However, 
increasing 7 above 10, will not change the diffusion con- 
stant much in the frequency range we discuss. 

The diffusion constants in ID and 3D media show the 
same behavior. Of course, the resonances at low velocity, 
are caused by the fact that all scatterers are assumed to 
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have equal size. When scatterer sizes (or velocities) are 
allowed to vary, the resonances will be averaged out. 

When the scatterer velocity is zero we obtain an im- 
penetrable model scatterer. This is not a useful model 
scatterer, as in the low frequency range the mean free 
path (and thus the diffusion constant) differ considerably 
from the penetrable scatterer case. The reason for this is 
that the limits for lu and 7 — > do not commute, as 



lim lim if = constant, 



while 



lim lim £f 

7^0 w^O 



00. 



(74) 



(75) 



The effect of this is that at the longer wavelengths, £f for 
the impenetrable scatterer is orders of magnitude smaller 
than if for non-zero values of 7. 



VII. CONCLUSIONS 

We have calculated the transport of energy and inten- 
sity in disordered ID, 2D and 3D (infinite) media emitted 
by a monochromatic source. Using the ladder approxi- 
mation to the Bethe-Salpeter equation we explicitly show 
that the total intensity is well approximated by the sum 
of the coherent and the fully developed diffuse wave field, 
for all source-receiver distances. Energy transport and 
intensity propagation in 2D disordered systems shows in- 
terestingly different behavior compared to the 3D case: 
In 2D, the average energy flux depends on the mean free 
path and the effective transport velocity depends differ- 
ently in terms of the scattering parameters. The (gradi- 
ent of the) intensity as a function of the source-receiver 
distance, on the other hand, behaves similarly in the 2D 
and the 3D case. The monochromatic source enables us 
to investigate the frequency dependence due to finite size 
scatterers on the macroscopic diffusion constant. For a 
monodisperse distribution of scatterers shape resonances 
show up in the relevant frequency range for low inter- 
nal scatterer velocities (7 small). In this frequency range 
(where scattering is expected to be isotropic) the depen- 
dence of the scattering properties on frequency can not 
be neglected. This means that descriptions of broadband 
pulse propagation through these media should in princi- 
ple incorporate both frequency dependent and multiple 
scattering effects. The development of a workable Ward 
identity in this case remains a challenge, however. Fi- 
nally, we want point out that our model describes trans- 
port of scalar acoustic waves but results can be extended 
and many conclusions should also apply to vector wave 
fields random media. 
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APPENDIX A: ENERGY AND INTENSITY IN 2D 

In this appendix we derive the configuration-averaged 
intensity and energy flux in a disordered 2D medium. 
Starting point is the 2D Green function propagator 

Q ^^ / -i^^o'^ ii^r + V {2if)) r) if ^ > 
I \Hl^\{^^r~i/{2if))r) if c^<0 ■ 



We use the properties 

ijf ((«.-z/(2^/))r) 



and 



= -z-Xo((W + l/(2^/))r), 



(Al) 



(A2) 



(A3) 



to express the Hankel functions {H^ ^) in terms of modi- 
fied Bessel function of the second kind {Kq). The Fourier 
transform of the coherent intensity 



no(fc;c^) = 27r/ |G(r;c^)r Jo(A;r), (A4) 



is then obtained from Ref. 24] and using properties of 
the associated Legendre polynomials 25] as 



Ho {k-Lo) 



with 



2U} 



(A5) 



(A6) 



This is can be rewritten as 



arcsin ( ^^^^pllz^^' 



TT 



^l + {kiff^{2K4ff^{klff 



(A7) 



no(fc;a;) is real, continuous and differentiable for all 
(real) fc > 0. 

The flux in the 2D system is given by 



(n • (r)) 



Q'^Pqoj C 
2 7' 



(A8) 
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where C is the constant to be calculated: 



drlm{G {r;uj)drG* {r;uj)}r^ 



C = 



The term in the denominator is easily obtained 



(A9) 



1 - 



{2Krif) 2Krif a.TCtan{2Kr£f) 



(AlO) 



The solution to the integral 

drlm{G{r;Lj)drG* {r;Lj)}r^ 



sgn (lu) 



drr^ lTji{{iKr + l/ {2£f)) 



xK„ {{-iKr + 1/ (2^/)) r) Ki {{iKr + 1/ (2^/)) r)} 



xF 2,2;3;l-ii±^^ 

\ {l-i2Krif) 



(1 - i2Kr£f) 



(All) 



can be found from Rcf. 24\. However, the proper solu- 
tion (on the right Riemann sheet) needs to be chosen in 
order to simplify the hypergeometric series F. One can 
check numerically that 



drlm{G{r;iu)drG* {r;Lu)}r^ 



sgn (lu) 

4^2£2 arctan(2K^£/) 



X 1 



1 



{2Kr(.f)^ 2Krif a.TCta.n{2Kr£f) 

Hence, C is given by 

sgn (w) 



(A12) 



C = 



47r2 



arctan (2Kr£f) . 



The intensity is proportional to the propagator 
n(r;ti;), expressed in terms of no(fc;w) by the Fourier 
transform (|46() . Only the gradient of the intensity is a 
well-defined property, so we calculate 



oo 

Vn(r;cc>) = -r 



fcVi (fcr) 



o'{k;Lo)-U^' {k = 0;u;) 



(A14) 

This part contains both the coherent and the scattered 
intensity. As the coherent intensity is known, we focus 
on the scattered intensity by calculating 



dk 

VH,, (r; w) = -f / — Ji (fcr) fc^n,, (fc; to) , (A15) 



with 



Use {k;uj) 



no-^(fc^O;c^)Ho {k;uj) 
n^'{k;u)-U^'ik = 0;u;y 



(A16) 



(|A15|I is the integral to calculate numerically when we 
want to calculate the gradient of the multiply scattered 
intensity. H^c {k; uj) is a monotonically decaying function 
with a maximum at fc = 0, that vanishes as /c — > oo. As 
the Bessel function is also decaying with r, we know that 
for r/£f > 1 



n,, (fc^O;cj) 
27rr 



(A17) 



and 



/ N ^arctan (2Kr^f) 1 i „ / . , „n 

WUtd (r) = -r I -g-^ {2nrlf) , (A18) 



where 



g{2K4f) = 1- 



{2Kr£f) 2Krif aTCtan{2Krif) 



(A19) 



(A13) td stands for "totally diffusive" . 
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